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We study fluctuation theorems for open quantum systems with a non-Markovian heat bath using 
the approach of quantum master equations and examine the physical quantities that appear in those 
fluctuation theorems. The approach of Markovian quantum master equations to the fluctuation 
theorems was developed by Esposito and Mukamel [Phys. Rev. E 73, 046129 (2006)]. We show that 
their discussion can be formally generalized to the case of a non-Markovian heat bath when the 
local system is linearly connected to a Gaussian heat bath with the spectrum distribution of the 
Drude form. We found by numerically simulating the spin-boson model in non-Markovian regime 
that the "detailed balance" condition is well satisfled except in a strongly non-equilibrium transient 
situation, and hence our generalization of the definition of the "entropy production" is almost 
always legitimate. Therefore, our generalization of the fluctuation theorem seems meaningful in 
wide regions. 
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I. INTRODUCTION 



The research of non-equihbrium statistical physics has been energeticahy developed for the last decade in the context 
of fluctuation theorems and the Jarzynski equality [H-Q- particular, the case of open quantum systems, where a 
reservoir with an infinite number of degrees of freedom is attached to a local quantum system with a finite number of 
degrees of freedom, is of great interest presently [1, l9l-[l6| . The most common approach to discussing the fluctuation 
theorems for an open quantum system is by the twice energy measurement on the local system together with the 
reservoir. There, one often assumes that the coupling between the local system and the reservoir is weak. This 
approach, since the measurements are done on both the local system and the reservoir, can be viewed as manipulating 
the total isolated system. 

In contrast, we here consider whether there exist any fluctuation theorems for an open quantum system that are 
expressed solely in terms of quantities of its local system. Fluctuation theorems of such kind were indeed proposed by 
Esposito and Mukamel under the context of Markovian quantum master equations. Their approach is a quantum 
analog to the fluctuation theorems of classical stochastic processes discussed by Crooks [13, 13 ^^'^ Seifert [l^. Note, 
however, that the entropy production that appears in their formalism is not a thermodynamic entropy production; it 
is an information entropy production. We will use the term entropy production in this sense in the following. 

In the present paper, we will first show that we can formally generalize the discussion by Esposito and Mukamel 
to the case of the dynamics in a non-Markovian heat bath. To this end, we employ a method called hierarchy 
equations of motion, which describes the reduced dynamics in a non-Markovian heat bath. The generalization to the 
non-Markovian dynamics in classical systems was already done in Refs. (20l - [2^ . 

When we discuss the fluctuation theorems, defining the work, the heat, and the entropy production is also an 
essential problem 0, [l^ [23l - l25| . As was done by Esposito and Mukamel [l^ > '^c can proceed with formal definitions 
based on quantum analogs to classical stochastic processes as working assumptions. However, it is necessary to examine 
whether the definitions are indeed appropriate quantum-mechanically. We will examine the validity of the "entropy 
production," the "entropy flow," and the "heat" that appear in the fluctuation theorems under our generalization 
by analyzing the spin-boson model, numerically in the non-Markovian regions and analytically in the limit where 
the Born-Markov approximation and the rotating-wave approximation are applicable. We will find that the "entropy 
production" indeed vanishes in equilibrium, and thus we can regard that our fluctuation theorem is legitimate even 
in non-Markovian regions except for strongly non-equilibrium transient states. 

As another point, if the "entropy flow" is related to the "heat" by the "microscopic reversibility," then we can 
rewrite our fluctuation theorem into the Crooks-type when the initial and the final states are the canonical states. 
We will find, however, that the "microscopic reversibility" is broken in wide parameter regions although it does hold 
in the limit where the Born-Markov approximation and the rotating-wave approximation are applicable in the driving 
protocols that we will investigate in the present paper. Furthermore, we will find that the "detailed balance" is almost 
satisfied in some parameter regions. 

The paper is organized as follows. In Sec. [Ill we will derive the fluctuation theorems for a quantum master equation 
in a non-Markovian heat bath. In this Section, we will assume that the heat bath is composed of the ensemble of 
harmonic oscillators that are linearly connected to the local system, but the Hamiltonian of the local system can 
be arbitrary. In Sec lIIIi we will show our numerical results for the spin-boson model; we will demonstrate the time 
evolution of the deviation from the "detailed balance." In equilibrium states, the "detailed balance" condition means 
that the "entropy production" is zero. We will also analyze the case where the local system is weakly driven by a 
Zeeman magnetic field periodically. The "entropy flow" and the "entropy production" are well-defined in this case, 
although they might be ill-defined when we drive the local system strongly or in different schedules. We will also 
examine whether any other useful relation holds; we will find that, although the "microscopic reversibility" which we 
will discuss in detail below, does not hold in general, there exist a case where the equation of the "detailed balance" 
instead is almost satisfied. 
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II. FLUCTUATION THEOREMS FOR A QUANTUM MASTER EQUATION IN A NON-MARKOVIAN 

HEAT BATH 



We consider a local system which is linearly connected to a Gaussian heat bath. The Hamiltonian of the total 
system reads 

H ^ Hs + Hb + H,^t, (1) 

i?B(&a,^D =I]^a&i&a, (2) 

a 

A„t = V^^Ca£a=:V^^c„(6, + 6t). (3) 

a a 

The Hamiltonian of the local system iJg is arbitrary. The heat-bath Hamiltonian H^i^a, h\) consists of the harmonic 
oscillators of an infinite number of modes Wq, where h]^ and ha are the creation and annihilation operators of each 
mode. We omitted the ground-state energy of the harmonic oscillators. The interaction between the system and the 
heat bath Hi-nt is linear in the position of the harmonic oscillator of each mode Xa = ^h/2maUJa + ba) with the 
coupling strength Cq, = Cq, \j2ma^al^- The operator V is an arbitrary one of the local system which depends on time 
in general. We approximate that the spectral distribution of the bath mode = '^^(c'^h/2maUJa)S{uj — cJq,) is of 
the Drude form: 

J{uj) = 7—^' (4) 

which is the Ohmic distribution with the Lorentzian cutoff 7 and a coefhcient C,. The cutoff 7 gives the decay rate 
of the canonical time correlation function of the heat bath, whereas C, is the coefficient related to the system-bath 
coupling strength Cq . For the region where the Born-Markov approximation is appropriate, i.e. , 



27r 

7: 



» C, (5) 



/3 

we refer to it as a Markovian heat bath; otherwise we refer to it as a non-Markovian heat bath. 

Esposito and Mukamel [l^ started their discussion with a Markovian quantum master equation of the form 

±\p{r)))^k{T)\p{T))), (6) 

where |p(t))) is the reduced density matrix of the local system and they used the following notation: 

\a,h))^\a){b\, (7) 
{{aMc.d)) = {a\c){d\b). (8) 

In ([HI), X is a dynamical semi-group acting on the reduced Liouville space of the local system. (We use the single- 
hatted letter as an operator acting on the Hilbert space and the double-hatted letter as an operator acting on the 
Liouville space.) The evolution of the reduced density matrix does not directly lead to the fluctuation theorem. 
In order to investigate the fluctuation theorem analogously to classical stochastic processes, they translated the 
Markovian quantum master equation in a form of a classical master equation using a time-dependent basis. They 
then constructed a "quantum trajectory" for the dynamics and discussed the forward and backward probabilities of 
the trajectory. The time-dependent basis is a basis which diagonalizes the reduced density matrix at each time, and 
hence the reduced density matrix is represented as 

(to^|/5(t)|77v) = {{m'^mr\p{T))) = PT{rn)5„i'm, 

l/5(T)))-^|m,)P.(m)(m,|, (9) 

m 

where we suppressed the subscript r for mr and on the right-hand side. Then we can regard the basis {I'ri.r)} as 
a set of states with probability {Pr{m)} at time r. By connecting these states, Esposito and Mukamel were able to 



construct a quantum trajectory (Fig. 1 in Ref. \XM) of the local system. Under this basis, they [IJ] had the classical 
master equation representation of the quantum master equation ([5]): 

y (Wr{m,m')Pr{m')-Wrim',m)Pr{m)), (10) 
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where 

Wr{m,m') = {(77iTmr\K{T)\m'^m/^)). (11) 

Wc will show that the discussion above by Esposito and Mukamel [3] can be formally generalized to the dynamics 
in the non-Markovian heat bath. In their discussion, it was essential that the equation of motion is expressed in a 
time-local form. In order to keep the time-local form for the non-Markovian heat bath, we employ a set of master 
equations called hierarchy equations of motion (26l - [29| . In the formalism of the hierarchy equations of motion, a state 
of the system is expressed with a set of infinite matrices in an extended space instead of a reduced density matrix. 
Thereby we can obtain a time-local equation of motion by taking account of the time correlation of the heat bath in 
the form of correlations among the matrices. 



A. Hierarchy equations of motion 



In the formalism of the hierarchy equations of motion |26l - l29j . a state of the local system is expressed as a set of 
an infinite number of matrices as follows: 

PoL.ir) ^ p^,%ir) ® pW...(r) ® • • • ^ pS-,,.,-i^) ® ' • ' 

pi%Ar) EE p(r), pi%,ir) ® p«,...(r) • ■ • {ai{r)}, (12) 

where /Sq^q (t) is the reduced density matrix p(r) in the usual sense, while {^i{t)} are the set of auxiliary matrices 
which possess the information of non-Markovian effects. These matrices follow from the equations of motion (|A15|) - 
(IA18[) in Appendix]^ The auxiliary matrices {^i{t)} are introduced for computational purposes only and possess 
no physical meaning themselves. The superscript n indicates the nth correction with respect to the decay rate of 
the heat bath 7 and the subscript {jk\k = 1,2,...} indicates the corrections with respect to the temperature, i.e. the 
Matsubara frequencies i^k = 27rfc//3fi.. We define its inner product as 

{{p{r);{Mr)}\p'ir);{alir)})) 

00 

EE Tr (pt(r)p'(r)) +^ Tr {aj{r)al{r)) . (13) 
1=1 

Then the hierarchy equations of motion (jA15|) - (jA18p can be formally expressed as 

^|p(r);{a,(r)})) =Le..(r)|/5(r);{a,(r)})), (14) 

where ihicr(''') is the generator that represents all the operations on the right-hand sides of Eqs. (jA15[) - (|A18[) . Note 
that Eq. ([Til) is derived under the assumption that the local system is linearly coupled to the Gaussian heat bath ([2]) 
with the spectrum distribution of the Drude form (jT]). This setting is often assumed when wc derive a Markovian 
quantum master equation microscopically. 



B. Birth-death master equation for the hierarchy equations of motion 

Analogously to Ref. [l^l, let {Iwr, riir] {6})} be the states that satisfy 

((m^,m^; {6}|p(r); {(7i(r)})) Tr \mr){mr\p{T)\mr){mr\ 

- Pr{m), (15) 

where {0} = {0, 0, . . . } denotes an infinite set of zero matrices in the space of the auxiliary matrices. Our purpose is 
to extract out the equation of motion for the probability Pr{m). Taking the inner product with {{rrirmr', {0}| in (|14p . 
we have 

((m.m.;{6}||-|p(r);{az(r)})) 

= ((m.m,; {6}|Zhicr(T)|/5(T); {a/(r)})). (16) 
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Next, wc decompose |/3(t); {a'i{T)})) into the following form: 



\piryAMr)}))=J2 



{{■mr\p{T)\m.r)) \mr){mr\ (X)]^(T((t) 



1=1 



|mrmr; {m{T)}))Pr{m). 



(17) 



In the first equality of (fT7|). we used the completeness relation for the reduced density matrix p{t). In the formalism 
of Ref. H, it reads 



\p{t))) ^^{\mr,nir)){{mr,mr\) |/5(t))) 
= ^ \mr){mr\p{T)\mr){mr\. 



Using (|17p in the right-hand side of (|T6| . we have 



{{m^mr;{6}\ — \p{T);{aiiT)})) 



The left-hand side of (fT9l) can be also recast as 



((m.r7v;{0}| — |p(r);{<7,(T)})) 



Tr (^|m,)(m.|Ap(r)^ + Tr (o ^^i{r)^ 
-^Pr{m) - ( ^("^^1"^^) ) Prim) 



dT 



Prim). 



Then we have 



Defining 



we can rewrite (1211) as 



dPr{m) 

dT 



= ^((m^TO^; {6}|£hicr(T)|m^m^; {(7;(T)}))P^(m'). 



W^(m, m';{(7/(r)}) = ((m^m^; {6}|£hicr(T)|TO^TO^; {(T;(r)})), 



— P,(m) = ^ VP^,(m, m'; {a/(T)}) P,(m'). 

Conservation of the probability during the evolution leads to 

0=|:Trp(r).A^P^(^) 



^ (Y,Wr{m,m'-{ai{T)})\ Pr{m'). 

m' \ m J 



Then we have 



^M^,(TO,m';{ai(T)}) = 0, 



(18) 



(19) 



(20) 



(21) 



(22) 



(23) 



(24) 



(25) 
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I.e. 



Wrim',m';{ai{T)}) 



m,m';{ai{T)}). 



(26) 



We thereby obtain an equation of motion for the hierarchy equations of motion in the form of a classical birth-death 
master equation: 



This is a generalization of Eq. ^TU\\ to the case of a non-Markovian heat bath. It is essential to consider an open 
system, since we can easily show that the "transition rate" Wr{Tn,m' ] {c^i{t)}) becomes zero for an isolated system; 
see Appendix [Bj Indeed, the transition of the local system occurs thanks to the coupling to the heat bath. 

Note, however, that the positivity of the "transition rate" WV(&T,aT) in (l22l) is not guaranteed. For the birth- 
death master equation in classical stochastic processes, we construct the equation of motion based on physically given 
transition rates; the transition rate there is positive by definition. In contrast, Eq. (|27p is the equation that we derived 
from the Liouville-von Neumann equation of the total system; Eq. (j27p does not necessarily have the meaning of the 
birth-death master equation in the sense of stochastic processes. 

Algebraically, Eq. (P7)) is correct in any parameter regions, except for the restriction that the local system is linearly 
coupled to the Gaussian heat bath with the spectrum distribution of the Drude form. Our calculation starts from 
the exact Liouville-von Neumann dynamics of the total system and uses no approximations before arriving at (j27p . 
Therefore, the simulations that we will show in Sec. IIIII are numerically exact. 



Equation (|27p naturally implies a quantum trajectory analogous to the trajectory in classical stochastic processes. 
We can regard the basis {|?TiT-)} as a set of states with probability {Pr^m)} and the local system hops from a state to 
another with the "transition rate" Wr{m' , m; {<Ji{t)}) at time r. We can construct a forward quantum trajectory by 
labeling the times r, (j = 1, 2, • • ■ , N) when the transitions occur and connecting these states (Fig. [T]): 



where rij represents the state after the jth transition at time Tj. Hence, ^(t) = nj for tj < t < Tj+i when j < N and 
n(T) = njv for tj^ < t < t. We set tq = and tn+i = t. 

Similarly, we can construct a backward trajectory corresponding to the forward one. The label of time for the 
backward trajectory is related to the forward one as tj O Tjv_j+i and the value (duration) of the time is related as 
Tj ~ t — TTv-j+i. The label of the state is related as = njv_j+i, so that the backward trajectory corresponding 
to n(r) is (Fig. [J) 




Wr{m\m;{ai{T)})Pr{m) . 



(27) 



C. Forward and Backward Trajectories 



n(r) = no — ?■ rii — )■ 77,2 • ■ • — ^ un, 



(28) 



^(t) = no — > ni — > 71,2 • ■ • — )■ 



^ TIN ^ nN^i nN-2 ■ ■ • no, 



(29) 



where we set tq = t and fj^+i = 0. 

The probabilities for the forward and backward quantum trajectories are obtained just by replacing WTj{nj,nj-i) 
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mJx) 













P(0) p(T,) p(tj p(T,) p(f) ^ 










r 


























«1 




► 



= T„ 



o=f, 



T3 ? = T4 



FIG. 1. (color online) A quantum trajectory n(r) for A'' = 3 in the case of a two-level system. Precisely speaking, the state 
mo and mi are time-dependent, and hence the axis should be dynamically modified, but it is suppressed here for simplicity. 
The labels with tilde represent the backward process. 



with Wr^{nj,nj^i; {ai{Tj)}) in Eqs. (62) and (66) of Ref. [ij: 



MF[n-(T)] = -Po(n-o) 

-N / 



dT'Wr'{m,n^-i;{aiiT')}) 



X Wr^{nj,nj-i;{ai{Tj)}) 



expl-^/ dT'Wr'(m,nN;{&iiT')})\ 

\ rn J 



MB ["-(?)] = ^o(?^o) 



N 



dT'Wr'{m,n^,i;{ai{T')}) 

1 

X Wr^{nj,nj.i;{ai{Tj)}) 



(30) 



(31) 



Hereafter, we set the probability of the final state of the forward trajectory equal to the probability of the initial state 
of the backward trajectory, i.e. Po("-o) = Pt{nN)- 

In the argument in Ref. [T3 | , the dynamics that gives the backward process was a hypothetical one; the real time- 
reversed dynamics does not satisfy the relation Wr(jh, fh') = Wt-rifn, fh') because the quantum master equation with 
the Markov approximation breaks the time reversal symmetry. In contrast, the backward process that we consider in 
(|3ip is truly the process of the time-reversed dynamics because the hierarchy equations of motion formally solve the 
total system and do not break the time reversal symmetry. Then we exactly have 



Wr{m,m';{ai{T)}) = m_?(m, to'; {(Ti(t - f )}). 



(32) 



Note that this condition dictates that the auxiliary matrices {o'lir)} also evolve backward because the dynamics of 
the heat bath must be also reversed. Therefore, precisely speaking, the backward process which satisfies also 
becomes a hypothetical one if we truncate the hierarchy, although it would be effectively a real time-reversed process 
if the conditions (|A11[) and (|A19p in Appendix |X] are satisfied. 
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D. Fluctuation theorems 



By repeating the same argument as in Ref. jl4| . we can show the fluctuation theorem for the dynamics in a 
non-Markovian heat bath. 

We define the entropy change As(t), the "entropy flow" Aso(t), and the "entropy production" Asi{t) along a single 
trajectory as 



As{t) 



Prin) 



■i(r) ] = 1 

= InPo(no) - InPt(nAr), 
Asc(t) = - / dr E 5(t - Tj) In 




--El- 



Asi{t) = 



Wr^{nj,nj_i;{ai{Tj)}) 
Wr^{nj-i,nj;{ai{Tj)})' 

Prin) 



V 
As(t) - Aso(t). 



(rtj ) Wr (nj - 1 , nj- ; {g; (r ) } ) 
P^(nj_i)Wr(7ij,nj„i; {o-;(r)}) 



(33) 



(34) 



(35) 



This type of definitions were first introduced by Schnakenberg [30| and used by Seifert [19| in order to discuss the 
fluctuation theorems in classical stochastic processes. Esposito and Mukamel [ij] extended these definitions for the 
classical dynamics to the ones for the Markovian quantum dynamics. Our definitions (|33P " P5p are their generalization. 

We also define the logarithm of the ratio of the probabilities for a forward trajectory and a backward trajectory as 



^J■B[n{T)\ 



(36) 



where F and B stand for the forward process and the backward process, respectively. Substituting (|30p and (|31l) into 
I, we have 



^^ HfHt)] ^ ^^ Po(no) _ y> ^^ Wr^{nj-i,nj;{ai{Tj)}) 
f^B[n{T)] PtiriN) j-[ Wr,{nj,nj^i;{ai{Tj)]) 

= As{t) - Asc(t) = Asi(t). 



(37) 



This relation immediately leads to the integrated fluctuation theorem 

n(T) 



1, 



(38) 



where the average (• • ■)f is the one over all forward trajectories. 

The detailed fluctuation theorem also holds. The probability that rF(i) is equal to a value of the "entropy produc- 
tion" n{t) is 



PYm)) = ((5(l)(t) - A.Si(t)))F 



(39) 



^(T) 
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then p7|) leads to our fluctuation theorem 

PFim) - E A^B[n(f)]e^^'(*)<5(17(i) - Asi{t)) 

n(r) 

= e^^«5]MB[«(T)]<5(n(0-rp(i)) 
= e^^(*)^MB[H(f)]<5(17(i) + rB(t)) 

n(r) 

=: e"(*)pB(-f^(i)), (40) 

where we used X]k(t) ~ En(T)- Although Eq. pO)) is formally the same as the fluctuation theorem for the case in a 
Markovian heat bath, the "transition rate" Wr{m, m'; {criir)}) contains the effect of non-Markovian properties in the 
auxiliary matrices {d'i{Tj)}. 

Note that the derivation of the fluctuation theorems here are totally formal; we need to examine whether the 
definition of the "entropy production" ([55]) is physically legitimate. 



E. Points to be checked 



It is crucial for a fluctuation theorem whether the quantities Asi(t) and Asc(t) that we deflncd above are indeed 
appropriate as the entropy production and the entropy flow, respectively. We here list two points to be checked in 

secnni 

First, as a property of the entropy production, we expect it to be zero in equilibrium. Therefore, Asi(t) defined in 
([55]) should vanish there, i.e. 

Pr{nj)Wr{nj.i,nf, {ai{T)}) ^ ^ 

which is nothing but a detailed balance condition. If this condition is satisfied and the "transition rate" Wr{m, m'; {o'i{t)}) 
is positive, then we can regard the formal equality (j37p as the fluctuation theorem. Otherwise, p7p is an equality 
without any meanings of entropies. As we mentioned at the end of Sec. IIIB| the "transition rate" VFr (m, m'; {cri{T)}) 
can be negative; we will numerically demonstrate in Sec. IIII Al that it indeed becomes negative during a strongly 
non-equilibrium transient state. In such a case, the dcflnition (j35|) of the entropy production as well as (j34p of the 
entropy flow are inappropriate because the inside of the logarithm of ((35|) and (p4|) . or the left-hand side of (|4ip . 
becomes negative. Remarkably, we will flnd in Sec. IIII Bl that except for the transient case, the "transition rate" 
Wr{m, m'; {o'z(t)}) is almost always positive and the "detailed balance" is almost satisfied. Furthermore, we will show 
in Appendix [Cl that the "transition rate" WV(m,m'; {^((t)}) is always positive and the equilibrium state satisfies the 
"detailed balance" in the limit where the Born-Markov approximation ([5]) and the rotating-wave approximation (|C1[) 
are applicable. 



Next, following Ref. [IJ], we define the "heat" that flows out of the local system along a trajectory by 

N 

Ut) ^ E {{n,\Hs{r,)\n,) - K_i|ffs(r,)|n,_i)) . (42) 

As a property of the entropy flow, we expect that it satisfles the following equality: 

Ase(i) = -/3 5s(i)- (43) 

After deflning q{nj, nj^i^Tj) ~ {rij \Hs{Tj)\nj) — {nj^i\Hs{Tj)\nj^i) to be qs{t) =: X^jLi $'(%'i '^j-ii T?)' '^^ recast 
our expectation (|43p to 

^' ^\\\ = exp [/3q{nj , _ i , r,- )] . (44) 

The equality (|44p is what we call the "microscopic reversibility." Note that this relation does not affect the validity of 
P7p . Although the probability distributions of the initial and the flnal states are arbitrary in general, when they are 
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both the canonical state, i.e. , Po{no) = exp[—(3{E{0) — F{0))] and Pt{npf) = cxp[—f3{E{t) — F{t))], we can rewrite 
([57| into the form of the Crooks-type fluctuation theorem: 



MB[n(r)] 

where AF = F{t) - F{0) and we defined the "work" 



''^["^^^l -eM~(3iwit)~AF)], (45) 



w{t) =AE~ q{t), (46) 

where AE = E{t) — E{0). Note that we can only check the consistency of the "entropy flow" Asc{t) and the "heat" 
q{nj, rij-i, tj); while the left-hand side of (|44p is not an established quantity as the ratio of the transition rates as we 
mentioned at the end of Sec. Ill Bl the "heat" q{nj , rij-i, Tj ) on the right-hand side is also a hypothetical quantity. In 
Sec. IIII Bl and in Appendix [Cl we will show that the equality (|44)) holds when the Born-Markov approximation ([5]) 
and the rotating-wave approximation (jCip are appropriate, although it does not hold in general. 

In the following section, we will numerically examine the "detailed balance" (|4T|) and the "microscopic reversibility" 

da. 



III. NUMERICAL EXAMINATION OF THE "DETAILED BALANCE 
AND THE "MICROSCOPIC REVERSIBILITY" 



In order to examine the above points numerically, let us consider the spin-boson model: 



-f i^int(V^,V'^&a,&^.) +i/cou„tor(^,V'^), (47) 

Hsil ^t) ^ = ^(^t^ _ ^^t), (48) 

HBiL,bi) = J2f^c.biba, (49) 

a 

HinS,i^\b^,bi) ^V{ipJ^)Y,(^c.{K +bi), (50) 

Vi^', 4'^) = Via^ + V^a, = ViiiJ;^ +4') + V^2(^^^ - H^). (51) 

(^,^^) = E '°r^'f' - (52) 

a ^ 



The operators cTx and CTz are the Pauli matrices and Vi and V2 are the corresponding coefficients. The local system 
Hamiltonian Hsii', "0^) is the two-level system with its energy difference fvujo- We added a counter-term i?counter('0, "0^) 
to the local system Hamiltonian in order to maintain the translation invariance of the system with respect to the heat 
bath [M [33. We set ?i = 1 below. 

In numerical simulation of the hierarchy equations of motion, the set of parameters {N, K) in (|Alip and (jA19[) (see 
Appendix[X| determines the accuracy of the calculation. As we take the values of {N, K) larger, the physical quantities 
converge to a certain value, which is numerically exact. The present paper is based on the program nonMarkovianOQ, 
which is distributed on the web site of Yoshitaka Tanimura |33| . 

We will mainly examine Eq. (|4ip in Sec. IIII Al and examine Eq. (|44p in Sec. IIII Bl while we will also investigate 
some other properties of the dynamics. We again stress that the calculation is numerically exact in the sense that we 
did not use any approximations in deriving Eq. (|14p. In the region where the Born-Markov approximation and the 
rotating wave approximation are applicable, we can analytically investigate the behavior with the quantum optical 
master equation (Appendix [Cjl . and therefore we will compare our results to those cases. 



A. Thermalization after the energy measurement and the "detailed balance" 

Let us consider the following protocol: we measure the energy of the isolated two-level system at r = and find 
it in the ground state. Then we connect the two-level system to the heat bath which is in the thermal equilibrium 
state at the inverse temperature /3 and let the total system evolve without any other perturbation. As long as the 



11 



coupling between the local two- level system and the bath contains a non- vanishing term of Vi, the local two- level 
system exchanges the energy with the heat bath and relaxes to a stationary state. Since we are not perturbing the 
total system with any external fields, we take this stationary state as the equilibrium state. Note that the equilibrium 
state of the two-level system differs from the canonical state in general; the total system should be in the canonical 
state in thermal equilibrium, but the reduced system is not. 

The relation that we examine here is the "detailed balance" (j4T|) . or equivalently 



Wria,b;{ai{T)}) ^ Prja) ^ (ar |jo(T)|a^) 

W^{b,a-{ai{T)}) PAb)- {br\pir)\br)' ^ ' 

The set of states {|ar}, |&t)} is the basis that diagonalizes the reduced density matrix at time t. We label the elements 
of the basis |ar) and \br) in the non-decreasing order of the eigenvalues of p{t). We denote the ratio of the "transition 
rate" between the states Or and br (the left-hand side of ([55)1 ) as 

. Wria,b;{ai{T)}) 

DBl(t) .= — — (54) 

Wr{b,a;{ai{T)}) 

while the ratio of the probability of the state at time t (the right-hand side of ([5^ 1 as 

= {bMr)\br)- ^''^ 

In general, the detailed balance condition is expected to hold in the case of Markovian dynamics in thermal equilibrium. 
We simulate the time evolution of the ratio DBl(t)/DBr(t) in order to examine whether such a condition holds in 
generic cases. We will show that the equality holds after the system reached the equilibrium state, although the 
"transition rate" ([H]) can be negative during the transient state. 

We examine the case of Cj^-coupling (i.e. Vi — 1 and V2 = 0) to the bath as well as the case of {(Tx + (t^ )-coupling 
(i.e. Vi ^ V2 ~ 1). Figures [5^ and[5)D show the time evolution of the ratio DBL(T)/DBr{ (t) for the parameter region 
where the the Born-Markov approximation ([5]) and the rotating-wave approximation (jCl[) arc relatively appropriate 
in the cases of (7a;-coupling and {crx + (T2)-coupling, respectively: 7 = IOOwq, 27r//3 ~ 120a;o, and C = O.OOlwo- 

For any cases of Ca;- coupling, the basis is time-independent because no off-diagonal element comes out if the state 
before the thermalization is diagonal. We can understand this from ()A3P in the Feynman- Vernon theory; because the 
coupling is linear in the bath coordinate and the bath is Gaussian, the transition always occurs in even orders of a^, 
and thus it causes no off-diagonal elements. In contrast, the case of the {ax + (T2)-coupling is the one where the basis 
can be time-dependent. The time evolution of DBl(t)/DBr(t) in Fig.[2l3 is fluctuating compared to Fig. [2^. This 
fluctuation is due to the time evolution of the basis. 

All the results of the simulations show that the "detailed balance" ([5^ holds once the system reaches the equilibrium 
state: r > 90 in Fig. [2^ and r > 9 in Fig. [2}3. The detailed balance is a sufficient condition, but not a necessary 
condition for the equilibrium state. Therefore, this is a highly nontrivial result. It means that the "entropy production" 
that we defined in p5|) satisfies the important requirement of the entropy production; it vanishes in the equilibrium 
state. 

Although everything seems fine in the equilibrium states, as we can see from the result of the simulation, 
Wr{rn,m';{ai{T)}) can be negative (i.e. DBL(T)/DBR (r) < 0) in transient states of the non-Markovian region. 
Hence we cannot always regard Wr{m,m') as the transition rate of the classical birth-death master equation. More 
importantly, it implies that the "entropy fiow" ([M)) and the "entropy production" ([55]) can be ill-defined because the 
arguments of the logarithms becomes negative. In contrast, we can prove that the "transition rate" Wrim^m!) is 
always positive in the case of the quantum optical master equation (Appendix [U]) , i.e. , the Born-Markov approxima- 
tion ^ and the rotating-wave approximation (|C1|) are applicable. For the results of Figs. [5^1. andl^b, we expect that 
([5|) and (|C1[) are appropriate as is shown above and hence the "transition rate" would be always positive. The system 
during the transient state, however, must be in a non-Markovian region after the strong perturbation of the energy 
measurement. This is indeed indicated by the fact that the ratio DBl(t)/DBf{,(t) changes drastically as we include 
more temperature corrections of the hierarchy equations of motion, i.e. , increasing the number of K. Except for such 
strongly non-equilibrium cases, our definitions of the entropy change, the entropy flow, and the entropy production 
seem legitimate as well as our fluctuation theorems (|37)) . ((38|) . and (|40)) . 



B. The "microscopic reversibility" and other possibihties of equahty 

In order to see whether the definition of the "entropy flow" p4p and the "heat" are legitimate, we examine the 
"microscopic reversibility" pij) . Just as the "detailed balance," the microscopic reversibility is a general relation in 
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FIG. 2. (color online) The evolution of the ratio DBL(T)/DBR(r) in the cases of (a) (Ji^-coupling and (b) {ax + (7z)-coupling 
after the energy measurement at r = 0. The parameter values are 7 — lOOojo, ficjo = 0.05, and ( = 0.001a;o, which is the case 
where we expect the Born-Markov approximation ^ and the rotating-wave approximation (IClf) to be relatively appropriate; 
min [7, 27r//3] ^ ^, and too 3> C- The time step is 10"'* in the simulation time. The lines indicate, from the top to the bottom, 
{N,K) = (4,1), (5,2), (6,3), (7,4), (8,5), (9,6), and (10,7). The line converges to the case where an infinite number of 
auxiliary matrices are taken into account, i.e. the full solution. 



the case of Markovian dynamics. We cannot expect such a relation to exist out of the Markovian limit; nevertheless, 
it is worth testing. We also investigate whether any quantities are related to the ratio of the "transition rate" DBl(t) 
other than the "microscopic reversibility" (|44| . 

The energy-measurement protocol of Scc. lIII Al was suitable to see the non-equilibrium behavior, but it is difficult 
to investigate the non-Markovian behavior because the convergence with respect to the hierarchy of the hierarchy 
equations of motion gets extremely slow. 

We here consider the following protocol instead: we prepare the total system in the thermal equilibrium state at 
time r = and drive the local two-level system by a sinusoidal external field for r > 0. In other words, we add to the 
Hamiltonian of the total system ([1]), the Hamiltonian of the external field 

H,^S, r) = h{T)a, = h{T){i^^ - U^) (56) 

with 

ft,(r) = Asin(a;cxtT), (57) 

where Wext is the angular frequency of the external field. Therefore, the local system does not go to equilibrium but 
to a stationary state at best. In the following, we consider the case of the {ax + (Tz )-coupling to the bath. 

In the present Subsection, we compare the ratio of the "transition rates" DBl(t) in ([51]) with the following three 
quantities. First, we compare it with DBR(r), the right-hand side of the "detailed balance" as we defined in (|55p . 
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The equality 

DBL(r) = DBR(r) (58) 



means that the "detailed balance" (|4ip is satisfied even in the non-equilibrium state. Note that one cannot expect 
the "detailed balance" to hold out of equilibrium in general. 
The second quantity that we compare with DBl(t) is 

UBr [ti[T)) .- ^^^|^eq(;,(^))|5;) ' l&y) 

where {|a^), \b'^)} is the basis that diagonalizcs the equilibrium reduced density matrix /5°'^(/i(t)). In actual calculation, 
we prepare the instantaneous equilibrium state /5°'i(/i(t)) by simulating the dynamics with /i(t) fixed. Although 

T)B^'^\h{T)) appears to be similar to DBr(t) in (|55p . we calculate the ratio of the equilibrium probability with the 
parameter ft.(r), the fixed value of the external parameter at time r. If the density matrix of the local system is equal 
to the equilibrium density matrix, we have 

DBR(r)lDB|,'='^)(/i(r)). (60) 

As is shown in Sec. lIII Al the "detailed balance" DBl(t) = DBr(t) should be satisfied in the equilibrium state, and 

thus we call that the system is in the quasi-equilibrium state when DBl(t) ~ DB^'^\h(T)) is satisfied. We can 

regard that the difference between DBl(t) and DB^^\h(T)) as an indicator of how close the system is to the thermal 
equilibrium. 

The third quantity that we compare with DBl(t) is MRr{ (T) := exp[— /3g(6, a, r)]. The equality 

DBl(t) =MRR,(r) (61) 

means the "microscopic reversibility" (|44p . As is shown in Appendix O the "microscopic reversibility" (|44p holds in 
the case of the quantum optical master equation, i.e. , when the Born-Markov approximation ([5]) and the rotating-wave 
approximation (jCip are applicable. 

Among these quantities that we compare, DBl(t), DBR(r). and DB^'^'' (/i(r)) are the quantities that are related to 
the birth-death master equation (j27p . while MRr(t) — exp[— /3q(6, a, r)] is a (hypothetical) thermodynamic quantity 
suggested by Esposito and Mukamel Q . 

In order to check the consistency, let us first consider the parameter region where the Born-Markov approximation 
([5]) and the rotating-wave approximation (jCip arc relatively appropriate: 7 = lOOwoi 27r//3 ~ 120wo, and C — O.OOIwq- 
The result in Fig.|3K shows that the state almost satisfies the "detailed balance" (|58| as well as ([60|) . so that the system 
is in a quasi-equilibrium state. It indicates that the break of the "microscopic reversibility" (|6T|) is very small, and 
thus the "heat" q{b,a,T) and the "entropy flow" Ase{t) should indeed satisfy the relation P5|) . or equivalently, are 
legitimate as the heat and the entropy flow in the limit of the Born-Markov approximation ([5]) and the rotating-wave 
approximation (jCip . The result is completely consistent with the result of the quantum optical master equation in 
Appendix O 

Figure is for the parameter region where the Born-Markov approximation ([5]) and the rotating-wave approxi- 
mation (jCip do not work at all: 7 = Swo, 27r//3 ~ 12^0, and C = wq. Although the frequency of the external field 
Woxt is not very fast compared to the typical relaxation time of the reduced density matrix, Fig.[3]3 indicates that 
the difference of DBp (/i(t)) from DBl(t) is of the same order as the amplitude of DB^'^-' (/i(t)) itself and thus the 
system is out of equilibrium significantly. 

It is then remarkable that the deviation of DBr(t) from DBl(t) is much smaller than that of DB^'^\h{T)) in the 
present case; even though the system is in a strongly non-equilibrium situation, not even in the quasi-equilibrium 
state, the "detailed balance" ([58]) is almost satisfied. Note that it does not mean the "entropy production" (|35|) is 
almost equal to zero since the system is out of equilibrium and ^-(n) in ([35]) is nonzero. The only case that we have 
found so far where the definition and ([55)1 is inappropriate in principle is the transient state in Sec. IIII Al 

Finally, in contrast to Fig. [3ji, MRr(t) largely differs from DBl(t), and thus the "microscopic reversibility" p4)) 
does not hold at all in our sense. It means that, in the region where the Born-Markov approximation ([Sj and the 
rotating-wave approximation (|Cip arc inappropriate, our definitions of the "heat" q{b, a, r) and the "entropy flow" 
Ase(t) are not legitimate and do not satisfy the Crooks- type fluctuation theorem 
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FIG. 3. (color online) (a) The time evolution of DBlCt) (solid line), DBR(r) (A dots), DBJ,"''' (/j(r)) (x dots), and MRR(r) 
(+ dots) for A = 0.25tJo and ujcxt ~ W~^. The time step is lO^'* in the simulation time. The horizontal axis is scaled as ojextT. 
The parameter values are 7 — lOOojo, Ptoo = 0.05, and = O.OIojq, which is also the case where the Born-Markov approximation 
and the rotating- wave approximation are appropriate; min [7, 2-k j fS\ ^ and cjo 2> C- (b) The time evolution of DBl(t) (solid 
line), DBR,(r) (dots), DB^''''' (/i(r)) (solid line), and MRR(r) (solid line) for A = 0.25ajo and ojoxt = 10"^. The time step is 10"" 
in the simulation time. The horizontal axis is scaled as ojoxtT. The parameter values are 7 = 5ajo, /3ajo = 0.5, and C, — ojq, which 
is the case where the Born-Markov approximation and the rotating-wave approximation are inappropriate. In both cases, the 
hierarchy of (A*', K) is well converged. 



IV. CONCLUSION 

We formally generalized the fluctuation theorems for the Markovian quantum master equation to the case of the 
master equation in a non-Markovian heat bath in the case of a system linearly coupled to a Gaussian heat bath 
with the spectrum distribution of the Drude form. These fluctuation theorems arc based on the analogy with the 
classical stochastic processes. Since it was unclear whether the complete analogy holds, we numerically investigated 
the properties of the birth-death master equation ((27)) . The numerical results of Sec. IIIII showed that the "detailed 
balance" (|4ip holds in equilibrium (Scc. lIIlX| . Therefore, the "entropy production" ([55)1 is indeed zero in equilibrium 
and our fluctuation theorems ([57)1 . psp . and PO)) are legitimate. Furthermore, remarkably, the "detailed balance" 
almost holds even in the weakly non-equilibrium state (Sec. IIIIBp . We do not have a complete interpretation about 
this fact yet. The only exception that we found is that the "transition rate" p2|) on the quantum trajectory can 
be negative in the strongly transient dynamics after the energy measurement, and thus we cannot regard it as the 
transition rate in the sense of the classical stochastic processes. For such a state, the "entropy flow" and the "entropy 
production" may be ill-dcfincd because the arguments of the logarithms of p4p and (|55|) become negative. 
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Prom the dynamics of the sinusoidally driven system by the Zeeman field, we confirmed that the "microscopic 
reversibihty" (|33|) does not hold in the parameter region out of the Born-Markov limit, although it does hold in the 
quantum optical master equation, which is the case where the Born-Markov approximation ([5]) and the rotating-wave 
approximation (jCl[) are applicable (Appendix [C|) . One of the possibilities that we can think of is that the definition 
of the "heat" (|42|) is inappropriate in the non-Markovian regions. 

In summary, the "entropy production" (|35p that appears in the fiuctuation theorem (j37[) is appropriate other than 
in the strongly transient state. We cannot, however, rewrite ([37]) into the form of the Crooks-type fiuctuation theorem 
with the definition of the "heat" (|^^ . A proper definition of the heat in the non-Markovian regions might exist so 
that the Crooks- type fluctuation theorem may hold. 

It is difficult to investigate a state very far from equilibrium numerically, because we require a large set of (N, K) 
in order to make the result converge. We need a more efficient methodology in order to explore the properties of the 
dynamics on the quantum trajectory beyond the parameter regions that we computed here. 
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Appendix A: Hierarchy equations of motion 

Although the hierarchy equations of motion are described well in Rcfs. [26l - l29j . we here give a brief review in order 
to make our description more self-contained and clarify the notations. 



1. The Feynman- Vernon theory 

The hierarchy equations of motion describe the dynamics of the reduced density matrix of an arbitrary local system 
which is linearly connected to a Gaussian heat bath, i.e. , the Hamiltonian of (Il])-(|3]). We assume that the initial state 
of the total density matrix is a product state and the density matrix of the heat bath is in the canonical state, i.e. , 

/5tot(0) = m ® ^e-^^^ Zb = Tree-^^-, (Al) 

where Tre denotes the trace with respect to the bath degrees of freedom. According to the Feynman- Vernon theory 
[s^. l35j|. the time evolution of an element of the reduced density matrix p{ip, V'l ''") '■— (V'|TrBPtot(''')|V') with the above 
initial state is expressed as 

ri'(r)=Tpi,Tp{r)=^t,^'{T)=ili'f,Tp'{T)=4'f _ 

X /_ _ _ _ DtpDt/jDt/yDip' 

J,/,(O)=i/'iXO)=i/'i.'0'(O)=i/J'i,)/"(O)=V,' 

X e^(^«['''''''l-^^l^'^'l) J-Fv[?,V',W,V'']p(?i,V'i',0), (A2) 

where ip and -0' are the coherent states, 6*3 is the action of the local system, and the functional is the influence 
functional which contains all the information of the heat bath. Using the fact that the heat bath is Gaussian and the 
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coupling to the local system is linear, the influence functional reads 



where 



exp-^ 



.2 , diO dTiV'ir,] 

" Jo JQ 



-i j dT2J{uj)sin[uj{Ti ~ T2)]V°{t2) 





+ / (iT2 J(w) cos[aj(Ti — T2)] coth 



X exp 



a 

yx(T) ■.^Vii',i;,T)-ViW,^',r), 



(A3) 



(A4) 

(A5) 
(A6) 



2. Equations of motion 

In the case of the spectrum distribution of the Drude form Q , the influence functional becomes 



J"fv[i/','0,V''i V'l = exp 



JQ 



dri / rfr2$(ri)9(T2)7e' 



-7(^1 -T2) 



r z*^ 

xTTexp / dn rfr2$(ri)*fe(T2)j^fee 
fc=i LJo Jo 

xcxp (^-'SlJ\r^V°{n)V^n)^ , 



(A7) 



$(t) = iF''(T), 



e(r) = 

*fe(T) 



^^V^°(r) + ^cot(^)^x(.) 



2 



(A8) 
(A9) 

(AlO) 



where = 2Trk//3h is the Matsubara frequency. In the form (jA7p of the influence functional, the time scale of each 
element is clear. Hence, for k {> K + 1) that satisfies vi^ ^ ujq, we may approximate 



J/fcC 



-Mr-r') ^ S{t-t') 



(k > K + 1). 



(All) 



In other words, the Markovian approximation is valid in the levels k > K + 1, but not in the lower levels. After some 
calculations, under the constraint that < /3h'y/2 < tt, we arrive at 

J"fv['0,'0,'0', V''] 



exp-i 



dn $(Ti)e-T^i (^-^ 'rfr2 7e(r2)eT^^^ +S(ri; 
n e^pj - I *(ri)e-'''=^i ^ ' dT2 iyk'^kiT2V''^'-^ + $(ri)^'fc(Ti 



(A12) 
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where we defined 



i-^cotr^\ 



V 2 y 



V- (r)FX (r) + ^^^^V°{r)V^ (r). 



(A13) 



Now we introduce the following matrices: 



(V'f,'0f,'r) = / dipidilJidi^'idipl 



■4>{'r)=ip!,i>{r)='4't,i>' {'r)=ip'!,4'' {t)=4!i 



-7T 



^(o)=V'i,'/'(o)=i/Ji,-i/j'(o)=-!/j'i,?/j'(o)='i/j; 
jk 



DiliDipDiji' Dip' 



dr' Uk-^kiT'V"^ 



(A14) 



The element with (n, ji, . . . , jx) = (Oj 0, . . . , 0) represents the original reduced density matrix. Hereafter, we omit the 
subscript f for the final state. We can express the derivative of the infiucncc functional in terms of these matrices, 
and thus the formal time derivative leads to the following hierarchy structure: 



d 



d 

^pS.!.,o(V','/''; 



9t 



K 



(A15) 



fc=l 



p^'.^..,o('/',^';^) 



(0) 



K 



E^'^o...,i,...,o(V',V'';r), 



(A16) 



fe=l 



Pi.o,...,o(V',V'V) = 



- l^zX + J^i + ^ *^fe + S j 



pi,o...,o('/':V'';T) 



if 



fe=l 

-i^l*lP0°.^..,0(^:V'V), 



(A17) 



if 



A;=l 



3l:---:3K 

(n) 



E'&p5?....,.+i,.....(^,^';r) 



fc=l 
if 



(n) 



(?/>,?/'';t)- 



(A18) 



fc=i 



In the above equations, $,0, ^'fe, and S are given in (|A8p . (|A9|) . (|A10p . and (|A13p with the replacement (t) — 
and where V"^ f = Vf - fV and = + fV . We can summarize ((M5|) - (|il8l) formally as 
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3. Terminators 



111 principle, the set of above equations contains the hierarchy that continues infinitely. Nevertheless, for the matrices 
with j2i ■ " ' Jk) that satisfy 



K 

nj + ^jki'k:>(^o, (A19) 



k=l 

we can approximate that such matrices obey the following equation of motion: 



dr 



k=\ 



This is written solely in terms of the matrix with (n, ji, . . . , ^md thus the hierarchy is truncated here. The 
matrices which obey (jA20[) are called the terminators. 
For the numerical calculation, we set 

K 

N^n + J23k» , (A21) 

^ mm (7,1/1) 

number that satisfies (jAlQp . 



4. Numerical implementation 

Analytical calculation of the hierarchy equations of motion turns out to be quite difficult in most cases, and thus 
it is usual to obtain the result with numerical simulations. In the above formulation, wc can obtain the exact result 
if we take N in (|A2ip and K in (jAlip infinite so that the approximations in (|A20p and (|A11[) may become exact. 
We cannot, however, implement this because it requires an infinite set of equations to solve. Hence we simulate with 
finite values of {N,K), which produces of course an approximate result. Nevertheless, if we increase N and K, the 
physical quantity that we are calculating converges to a certain value. We can regard this value as the numerically 
exact result since the simulations with the higher set of {N, K) would only give negligible differences. 



Appendix B: Isolated system 

We will show in the present Appendix that the "transition rate" Wr{m,m') vanishes for isolated systems. In the 
case of an isolated quantum system, we do not need to use the method of the hierarchy equations of motion, and 
thus we take over the formalism by Esposito and Mukamel [l3| (see Eqs. (|5))-([TT|) in Sec. [ll]). The system obeys the 
Liouville-von Neumann equation 

^\p{r))) = l\pir))} = -^[Hir),pir)] (Bl) 

and the corresponding "transition rate" is 

Wrim, m') 
= {{mr\C\m'^)) 



—iy{mr\H\7Ji'^){'m'^\mr) — {{mr\H\m'^){m'^\'mr)) 
0. (B2) 



Thus, no transitions take place. The picture of the quantum trajectory appears obviously due to the interaction with 
the reservoir. 
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Appendix C: "Detailed balance" and "Microscopic reversibility" 
of the quantum optical master equation 

In order to compare with our results in Sec. Illli we will show that the "detailed balance" (|41]) . or equivalently 
and show that the "microscopic reversibility" (pi)) holds for the spin-boson model within the Born-Markov 
approximation ([5]) and the rotating-wave approximation. The condition of the rotating-wavc approximation reads 

^0 » C, (CI) 

where ujq is the energy spacing of the local system and C is the coefficient, which is related to the system-bath coupling 
strength Cq in Eq. ([3]). These approximations are common in studies of the reservoir of the photon field, and hence 
we here treat such a case. As we did in Sec. IIIIl we discuss in the Schrodingcr picture and set H = 1. 



1. Quantum optical master equation and the ratio of the "transition rate" 

For the Hamiltonian of the local system of Hg ~ uj{t^3/2, the equation of motion for the reduced density matrix 
p{t) is given by the quantum optical master equation [32| 

+ ld{T){N + 1) ^o-_/5(r)<T+ - i(T+(7_/5(r) - i/5(r)o-+o-_ 
+ 7d(r)iV (^&+p{r)<J- - ^a-a+p{T) - ^p{T)a^a+^ (C2) 

= : -i[Hi^s,p{r)]+i>\pir))) iLs+ t>Mr)))- (C3) 

In deriving (|C3[) . we used the Born-Markov approximation ([5|) and the rotating-wave approximation (jCip . The 

operator Hls is the Hamiltonian with the Lamb-Stark shift and £ls is its Liouvillian, 2? is referred to as the dissipator, 
(7+ = ((7a; -I- icTy)/2, a- = [dx - io-y)/2, and 

Ur) = Mli;^, N := iV(.(.)) = ,]^. . (C4) 

3c'^ exp(/3cj(T)j - 1 

The factor c is the speed of light and d is the off-diagonal clement of the dipolc operator, 

^dipoic(T) = d(T_e-'"(^)- + fTa+e'-^t^^^ (C5) 

This corresponds to the case of cr^-coupling in Sec. IIIIl In the current model and approximations, the (ctj; +(72 )-coupling 
would only cause the energy shift from the case of (Ta;- coupling. 

Since Eq. (|C3|) is of the form ©, we follow Esposito and Mukamel [13] here. The general solution of the density 
matrix is expressed as 

p(r) = e-'^-^^^i(l + (a(T)) ■ a)c""^^ 

^ ^-^rH.S A(l + (^3(^))) , ^ /M^)) ^ c^^- . (C6) 

The angular bracket means (■•■)= Trs(- ■ • /5(t)), where Trs is the trace over the local system. Each of these averaged 
values evolves as 

d 7d(T)[27V + l] 

^(cri(T)) = (cri(r)), (C7) 

d 7d(T)[2jV + l] 

^('^3(r)) = -7d(r)[27V + l](a3(r)} -7d(r). (C9) 
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Since the coupling is weak under the present approximations, the initial equilibrium state is well approximated as 



4/3(^(0) Q 



Thus, the density matrix in the interaction picture is a diagonal matrix all the time, which means that the time- 
dependent basis (jlSp becomes time independent: 

K>=(S) = (i)' = 

Because {{brlC-LsWr)) = as we see from (|B2[) . the ratio of the "transition rate" Wr{hT,aT) in ([TT]) reads 

Wr{hr,ar) = {{hr\V\ar)) (C12) 

= 7diV(|ai|2|62p + |a2nfeiH +7d|a2nfeir (C13) 

The second term is asymmetric with respect to the process reversal ar hr- We also see that the "transition rate" 
Wribrydr) is real and positive unlike the general expression ((221) . Substituting the specific form N = (e'^'^(^) - 1)-^ 
into ([UT3)) . we have for Eq. dM]) 

DB,(r) = ^ 
^ ' Wrib,a) 

e^-M|ai|2|62|2 + |a2p|6i|2 ' ^ ' 



2. "Detailed balance'' 



Let us consider the time evolution from the ground state of the two-level system with uj{t) = wq, just as in Sec. lIII K\ 
According to (jC6[) " (|C9[) . the time evolution of the density matrix /5(r) reads 



4(1 - 



with 



Thus we have for Eq. ([54 



2N{ujo) + 1 2N{ujo) + 1 



T^p _ Mp{T)\ar) _ 1 + {(Tsir)) , 
UBR(,rj — — -. (C17) 

{br\p{T)\br) 1-((73(t)) 

As T — !• cx), we have DBr (t) exp{— /3ujo) . Therefore, the "detailed balance" (HJ), or ([55)) holds in equilibrium under 
the present approximations. 

3. "Microscopic reversibility" 

For the case where the local system is driven by the Zeeman magnetic field as in Sec. IIIIBl the right-hand side of 
([ii[) reads 

{br\Hs{r)\br) = {ar\Hs{r)\ar) = (C18) 

exp{/3[{br\Hs{T)\br) - {ar\Hs{T)\ar)]) = c-^^^^^\ (C19) 
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Comparing (|C19|) with (jC14p . we see that the "microscopic reversibihty" (|44l) holds in this case. Note, however, that 
it does not generally hold even if we use the Born-Markov approximation and the rotating- wave approximation, when 
we drive the system with the external field other than the Zeeman magnetic field. 
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